In [1]:
from sympy import *
from sympy.vector import CoordSys3D
N = CoordSys3D('N')
x1, x2 = symbols("x_1 x_2")
alpha1, alpha2 = symbols("alpha_1 alpha_2")
R, L, ga, gv = symbols("R L g_a g_v")
init_printing()
In [2]:
a1 = pi / 2 + (L / 2 - alpha1)/R
x = R * cos(a1)
y = R * sin(a1)
r = x*N.i + y*N.j
In [3]:
r
Out[3]:
In [4]:
v_temp = r.diff(alpha1)
dr_len = v_temp.magnitude()
v = v_temp / dr_len
v = trigsimp(v)
v_temp
Out[4]:
In [5]:
trigsimp(dr_len)
Out[5]:
In [6]:
n_temp = v.diff(alpha1)
k=trigsimp(n_temp.magnitude())
n = n_temp/k
q=1/(R*sqrt(1/R**2))
n = trigsimp(n).subs(q, 1)
n
Out[6]:
In [7]:
v.dot(n)
Out[7]:
In [8]:
n.dot(v)
Out[8]:
In [9]:
sympify(k)
Out[9]:
In [10]:
n.diff(alpha1)
Out[10]:
$ \frac { d\vec{n} } { d\alpha_1} = -\frac {1}{R} \vec{v} = -k \vec{v} $
In [11]:
v.diff(alpha1)
Out[11]:
$ \frac { d\vec{v} } { d\alpha_1} = \frac {1}{R} \vec{n} = k \vec{n} $
$ \vec{u} = u_v \vec{v} + u_n\vec{n} $
$ \frac { d\vec{u} } { d\alpha_1} = \frac { d(u_v\vec{v}) } { d\alpha_1} + \frac { d(u_n\vec{n}) } { d\alpha_1} = \frac { du_n } { d\alpha_1} \vec{n} + u_n \frac { d\vec{n} } { d\alpha_1} + \frac { du_v } { d\alpha_1} \vec{v} + u_v \frac { d\vec{v} } { d\alpha_1} = \frac { du_n } { d\alpha_1} \vec{n} - u_n k \vec{v} + \frac { du_v } { d\alpha_1} \vec{v} + u_v k \vec{n}$
Then $ \frac { d\vec{u} } { d\alpha_1} = \left( \frac { du_v } { d\alpha_1} - u_n k \right) \vec{v} + \left( \frac { du_n } { d\alpha_1} + u_v k \right) \vec{n}$
$ \frac { d\vec{u} } { d\alpha_2} = \frac { d(u_n\vec{n}) } { d\alpha_2} + \frac { d(u_v\vec{v}) } { d\alpha_2} = \frac { du_n } { d\alpha_2} \vec{n} + u_n \frac { d\vec{n} } { d\alpha_2} + \frac { du_v } { d\alpha_2} \vec{v} + u_v \frac { d\vec{v} } { d\alpha_2} = \frac { du_n } { d\alpha_2} \vec{n} + \frac { du_v } { d\alpha_2} \vec{v} $
In [12]:
R_alpha=r+alpha2*n
R_alpha
Out[12]:
In [13]:
R1=R_alpha.diff(alpha1)
R2=R_alpha.diff(alpha2)
In [14]:
trigsimp(R1)
Out[14]:
In [15]:
R2
Out[15]:
$ A = \left( \begin{array}{cc} \frac{\partial x_1}{\partial \alpha_1} & \frac{\partial x_1}{\partial \alpha_2} \\ \frac{\partial x_2}{\partial \alpha_1} & \frac{\partial x_2}{\partial \alpha_2} \end{array} \right)$
$ \left[ \begin{array}{cc} \vec{R}_1 & \vec{R}_2 \end{array} \right] = \left[ \begin{array}{cc} \vec{e}_1 & \vec{e}_2 \end{array} \right] \cdot \left( \begin{array}{cc} \frac{\partial x_1}{\partial \alpha_1} & \frac{\partial x_1}{\partial \alpha_2} \\ \frac{\partial x_2}{\partial \alpha_1} & \frac{\partial x_2}{\partial \alpha_2} \end{array} \right) = \left[ \begin{array}{cc} \vec{e}_1 & \vec{e}_2 \end{array} \right] \cdot A$
$ \left[ \begin{array}{cc} \vec{e}_1 & \vec{e}_2 \end{array} \right] = \left[ \begin{array}{cc} \vec{R}_1 & \vec{R}_2 \end{array} \right] \cdot A^{-1}$
In [16]:
m11=R1.dot(N.i)
m12=R2.dot(N.i)
m21=R1.dot(N.j)
m22=R2.dot(N.j)
A=Matrix([[m11, m12], [m21, m22]])
A
Out[16]:
In [17]:
A_inv = trigsimp(A**-1)
sympify(trigsimp(Matrix([R1, R2]).T*A_inv))
Out[17]:
In [18]:
trigsimp(A.det())
Out[18]:
In [19]:
g11=R1.dot(R1)
g12=R1.dot(R2)
g21=R2.dot(R1)
g22=R2.dot(R2)
G=Matrix([[g11, g12],[g21, g22]])
G=trigsimp(G)
G
Out[19]:
In [20]:
G_inv = G**-1
In [21]:
dR1dalpha1 = trigsimp(R1.diff(alpha1))
dR1dalpha1
Out[21]:
$ \frac { d\vec{R_1} } { d\alpha_1} = \frac {1}{R} \left( 1-\frac{\alpha_2}{R} \right) \vec{R_2} $
In [22]:
dR1dalpha2 = trigsimp(R1.diff(alpha2))
dR1dalpha2
Out[22]:
$ \frac { d\vec{R_1} } { d\alpha_2} = -\frac {1}{R} \frac {1}{1-\frac{\alpha_2}{R}} \vec{R_1} $
In [23]:
dR2dalpha1 = trigsimp(R2.diff(alpha1))
dR2dalpha1
Out[23]:
$ \frac { d\vec{R_2} } { d\alpha_1} = -\frac {1}{R} \frac {1}{1-\frac{\alpha_2}{R}} \vec{R_1} $
In [24]:
dR2dalpha2 = trigsimp(R2.diff(alpha2))
dR2dalpha2
Out[24]:
$ \frac { d\vec{R_2} } { d\alpha_2} = \vec{0} $
$ \vec{u} = u^1 \vec{R_1} + u^2\vec{R_2} $
$ \frac { d\vec{u} } { d\alpha_1} = \frac { d(u^1\vec{R_1}) } { d\alpha_1} + \frac { d(u^2\vec{R_2}) } { d\alpha_1} = \frac { du^1 } { d\alpha_1} \vec{R_1} + u^1 \frac { d\vec{R_1} } { d\alpha_1} + \frac { du^2 } { d\alpha_1} \vec{R_2} + u^2 \frac { d\vec{R_2} } { d\alpha_1} = \frac { du^1 } { d\alpha_1} \vec{R_1} + u^1 \frac {1}{R} \left( 1-\frac{\alpha_2}{R} \right) \vec{R_2} + \frac { du^2 } { d\alpha_1} \vec{R_2} - u^2 \frac {1}{R} \frac {1}{1-\frac{\alpha_2}{R}} \vec{R_1}$
Then $ \frac { d\vec{u} } { d\alpha_1} = \left( \frac { du^1 } { d\alpha_1} - u^2 \frac {1}{R} \frac {1}{1-\frac{\alpha_2}{R}} \right) \vec{R_1} + \left( \frac { du^2 } { d\alpha_1} + u^1 \frac {1}{R} \left( 1-\frac{\alpha_2}{R} \right) \right) \vec{R_2}$
$ \frac { d\vec{u} } { d\alpha_2} = \frac { d(u^1\vec{R_1}) } { d\alpha_2} + \frac { d(u^2\vec{R_2}) } { d\alpha_2} = \frac { du^1 } { d\alpha_2} \vec{R_1} + u^1 \frac { d\vec{R_1} } { d\alpha_2} + \frac { du^2 } { d\alpha_2} \vec{R_2} + u^2 \frac { d\vec{R_2} } { d\alpha_2} = \frac { du^1 } { d\alpha_2} \vec{R_1} - u^1 \frac {1}{R} \frac {1}{1-\frac{\alpha_2}{R}} \vec{R_1} + \frac { du^2 } { d\alpha_2} \vec{R_2} $
Then $ \frac { d\vec{u} } { d\alpha_2} = \left( \frac { du^1 } { d\alpha_2} - u^1 \frac {1}{R} \frac {1}{1-\frac{\alpha_2}{R}} \right) \vec{R_1} + \frac { du^2 } { d\alpha_2} \vec{R_2}$
$\nabla_1 u^1 = \frac { \partial u^1 } { \partial \alpha_1} - u^2 \frac {1}{R} \frac {1}{1-\frac{\alpha_2}{R}}$
$\nabla_1 u^2 = \frac { \partial u^2 } { \partial \alpha_1} + u^1 \frac {1}{R} \left( 1-\frac{\alpha_2}{R} \right) $
$\nabla_2 u^1 = \frac { \partial u^1 } { \partial \alpha_2} - u^1 \frac {1}{R} \frac {1}{1-\frac{\alpha_2}{R}}$
$\nabla_2 u^2 = \frac { \partial u^2 } { \partial \alpha_2}$
$ \nabla \vec{u} = \left( \begin{array}{cc} \nabla_1 u^1 & \nabla_1 u^2 \\ \nabla_2 u^1 & \nabla_2 u^2 \end{array} \right)$
In [25]:
u1=Function('u^1')
u2=Function('u^2')
u1_nabla1 = u1(alpha1, alpha2).diff(alpha1) - u2(alpha1, alpha2) / R * (S(1)/(1-alpha2/R))
u2_nabla1 = u2(alpha1, alpha2).diff(alpha1) + u1(alpha1, alpha2) / R * ( 1-alpha2/R)
u1_nabla2 = u1(alpha1, alpha2).diff(alpha2) - u1(alpha1, alpha2) / R * (S(1)/(1-alpha2/R))
u2_nabla2 = u2(alpha1, alpha2).diff(alpha2)
# $\nabla_2 u^2 = \frac { \partial u^2 } { \partial \alpha_2}$
grad_u = Matrix([[u1_nabla1, u2_nabla1],[u1_nabla2, u2_nabla2]])
grad_u
Out[25]:
In [26]:
q=Symbol('q')
grad_u_down=grad_u.subs(1-alpha2/R, q)*G.subs((R-alpha2)/R,q)
#grad_u_down=grad_u*G
expand(simplify(grad_u_down))#.subs((R-alpha2)/R, q)
Out[26]:
$ \left( \begin{array}{c} \nabla_1 u_1 \\ \nabla_2 u_1 \\ \nabla_1 u_2 \\ \nabla_2 u_2 \end{array}
\left( \begin{array}{c} \left( 1-\frac{\alpha_2}{R} \right)^2 \frac { \partial u^1 } { \partial \alpha_1} - u^2 \frac {\left( 1-\frac{\alpha_2}{R} \right)}{R} \\ \left( 1-\frac{\alpha_2}{R} \right)^2 \frac { \partial u^1 } { \partial \alpha_2} - u^1 \frac {\left( 1-\frac{\alpha_2}{R} \right)}{R} \\ \frac { \partial u^2 } { \partial \alpha_1} + u^1 \frac {\left( 1-\frac{\alpha_2}{R} \right)}{R} \\ \frac { \partial u^2 } { \partial \alpha_2} \end{array} \right) $
$ \left( \begin{array}{c} \nabla_1 u_1 \\ \nabla_2 u_1 \\ \nabla_1 u_2 \\ \nabla_2 u_2 \end{array}
\left( \begin{array}{cccccc} 0 & \left( 1-\frac{\alpha_2}{R} \right)^2 & 0 & -\frac {\left( 1-\frac{\alpha_2}{R} \right)}{R} & 0 & 0 \\ -\frac {\left( 1-\frac{\alpha_2}{R} \right)}{R} & 0 & \left( 1-\frac{\alpha_2}{R} \right)^2 & 0 & 0 & 0 \\ \frac {\left( 1-\frac{\alpha_2}{R} \right)}{R} & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ \end{array} \right) \left( \begin{array}{c} u^1 \\ \frac { \partial u^1 } { \partial \alpha_1} \\ \frac { \partial u^1 } { \partial \alpha_2} \\ u^2 \\ \frac { \partial u^2 } { \partial \alpha_1} \\ \frac { \partial u^2 } { \partial \alpha_2} \\ \end{array} \right) $
In [51]:
from sympy import MutableDenseNDimArray
C_x = MutableDenseNDimArray.zeros(3, 3, 3, 3)
for i in range(3):
for j in range(3):
for k in range(3):
for l in range(3):
elem_index = 'C^{{{}{}{}{}}}'.format(i+1, j+1, k+1, l+1)
el = Symbol(elem_index)
C_x[i,j,k,l] = el
C_x
Out[51]:
In [53]:
C_x_1 = MutableDenseNDimArray.zeros(3, 3, 3, 3)
def getCIndecies(index):
if (index == 0):
return 0, 0
elif (index == 1):
return 1, 1
elif (index == 2):
return 2, 2
elif (index == 3):
return 0, 1
elif (index == 4):
return 0, 2
elif (index == 5):
return 1, 2
for i in range(3):
for j in range(i, 3):
for k in range(3):
for l in range(k, 3):
elem_index = 'C^{{{}{}{}{}}}'.format(i+1, j+1, k+1, l+1)
el = Symbol(elem_index)
C_x_1[i,j,k,l] = el
C_x_1[i,j,l,k] = el
C_x_1[j,i,k,l] = el
C_x_1[j,i,l,k] = el
if (i >= k or j >= l):
C_x_1[k,l,i,j] = el
C_x_1[k,l,j,i] = el
C_x_1[l,k,i,j] = el
C_x_1[l,k,j,i] = el
C_x_1
Out[53]:
In [ ]:
# for i in range(3):
# for j in range(i, 3):
# for k in range(3):
# for l in range(k, 3):
# elem_index = 'C^{{{}{}{}{}}}'.format(i+1, j+1, k+1, l+1)
# el = Symbol(elem_index)
# C_x[i,j,k,l] = el
# C_x[i,j,l,k] = el
# C_x[j,i,k,l] = el
# C_x[j,i,l,k] = el
# C_x
# for i in range(3):
# for j in range(3):
# for k in range(3):
# for l in range(3):
# elem_index = 'C^{{{}{}{}{}}}'.format(i+1, j+1, k+1, l+1)
# el = Symbol(elem_index)
# C_x[i,j,k,l] = el
# C_x[k,l,i,j] = el
# if (s < 3 and t < 3):
# elif (s==t):
# i,j = getCIndecies(s)
# k,l = getCIndecies(t)
# elem_index = 'C^{{{}{}{}{}}}'.format(i+1, j+1, k+1, l+1)
# el = Symbol(elem_index)
# C_x[i,j,k,l] = el
# def getCIndecies(index):
# if (index == 0):
# return 0, 0
# elif (index == 1):
# return 1, 1
# elif (index == 2):
# return 2, 2
# elif (index == 3):
# return 0, 1
# elif (index == 4):
# return 0, 2
# elif (index == 5):
# return 1, 2
# def getCalpha(C, A, q, p, s, t):
# res = S(0)
# for i in range(3):
# for j in range(3):
# for k in range(3):
# for l in range(3):
# res += C[i,j,k,l]*A[q,i]
# return simplify(trigsimp(res))
# C_alpha = MutableDenseNDimArray.zeros(3, 3, 3, 3)
# C_alpha_empty = MutableDenseNDimArray.zeros(3, 3, 3, 3)
# m11=R1.dot(N.i)
# m12=R2.dot(N.i)
# m21=R1.dot(N.j)
# m22=R2.dot(N.j)
# A=Matrix([[m11, m12, 0], [m21, m22, 0], [0,0,1]])
# A_inv=A**-1
# for s in range(6):
# for t in range(s, 6):
# if (s < 3 and t < 3):
# i,j = getCIndecies(s)
# k,l = getCIndecies(t)
# elem_index = 'C^{}{}{}{}'.format(i+1, j+1, k+1, l+1)
# el = Symbol(elem_index)
# C_x[i,j,k,l] = el
# C_x[k,l,i,j] = el
# elif (s==t):
# i,j = getCIndecies(s)
# k,l = getCIndecies(t)
# elem_index = 'C^{}{}{}{}'.format(i+1, j+1, k+1, l+1)
# el = Symbol(elem_index)
# C_x[i,j,k,l] = el
# for i in range(3):
# for j in range(3):
# for k in range(3):
# for l in range(3):
# c = getCalpha(C_x, A_inv, i, j, k, l)
# C_alpha[i,j,k,l] = c
# if (c != 0):
# elem_index = 'C_{{\alpha}}^{}{}{}{}'.format(i+1, j+1, k+1, l+1)
# el = Symbol(elem_index)
# C_alpha_empty[i,j,k,l] = el
# # trigsimp(C_alpha[0,2,0,2])
# C_alpha_empty
In [31]:
C_x
Out[31]:
In [ ]:
def contraction3DSameRank(A,B):
res = S(0)
for i in range(3):
for j in range(3):
res += A[i,j]*B[j,i]
return res
def contraction3D(C,e):
res = MutableDenseNDimArray.zeros(3, 3)
for i in range(3):
for j in range(3):
res[i,j] = S(0)
for k in range(3):
for l in range(3):
res[i,j] += C[i,j,k,l]*e[k,l]
return res
In [ ]:
e11 = Symbol("e_{11}")
e12 = Symbol("e_{12}")
e22 = Symbol("e_{22}")
e13 = Symbol("e_{13}")
e23 = Symbol("e_{23}")
e33 = Symbol("e_{33}")
# s11 = Symbol("s_{11}")
# s12 = Symbol("s_{12}")
# s22 = Symbol("s_{22}")
e=Matrix([[e11, e12, e13], [e12, e22, e23], [e13, e23, e33]])
s=contraction3D(C_alpha, e)
E=contraction3DSameRank(s, e)
In [ ]:
# e_alpha=G*A_inv*e*A_inv.T*G
# #e_alpha=A*e*A.T
# s_alpha=A_inv*s*A_inv.T
# E_alpha=contraction2D(s_alpha, e_alpha)